A novel frequency-dependent lattice Boltzmann model with a single force term for electromagnetic wave propagation in dispersive media

Electromagnetic wave simulation is of pivotal importance in the design and implementation of photonic nano-structures. In this study, we developed a lattice Boltzmann model with a single extended force term (LBM-SEF) to simulate the propagation of electromagnetic waves in dispersive media. By reconstructing the solution of the macroscopic Maxwell equations using the lattice Boltzmann equation, the final form only involves an equilibrium term and a non-equilibrium force term. The two terms are evaluated using the macroscopic electromagnetic variables and the dispersive effect, respectively. The LBM-SEF scheme is capable of directly tracking the evolution of macroscopic electromagnetic variables, leading to lower virtual memory requirement and facilitating the implementation of physical boundary conditions. The mathematical consistency of the LBM-SEF with the Maxwell equations was validated by using the Champman-Enskog expansion; while three practical models were used to benchmark the numerical accuracy, stability, and flexibility of the proposed method.


Results
To validate the proposed LBM-SEF model, we have carried out three typical cases which have either analytical or numerical exact solutions. The simulation results indicate that the new model is capable of accurately reproducing the EM wave propagation in both dispersive and non-dispersive media. Subsequently, we proceed to create more intricate simulations.

EM pulses in a non-dispersive media.
To conduct the initial benchmark, we conduct a simulation of a Gaussian pulse traveling through a one-dimensional array of L cells with periodic boundary conditions and crossing a dielectric interface. The first half of the simulation space, x < L/2 , represents a vacuum with a permittivity of ε 0 , while the second half, x > L/2 , represents a non-dispersive medium with a relative permittivity constant of ε r = 2.0 . The initial pulse is characterized according to previous studies 12,30 : where α is the pulse width, E M and H M are the amplitudes of electric and magnetic field of the pulse, respectively. We choose L = 800 , E M = 1000 , α = 30, and x c = 250 . And continues boundaries are applied at the boundaries. The initial condition and the electric field at t = 300 are shown in Fig. 1 Based on the work of Dhuri et al. 32 , we delved deeper into the dispersion properties of the 1D LBM-SEF model by examining the response to plane waves in homogeneous, non-dispersive media. Figure 2 compares the dispersion relations generated by our model with those of previous FDTD 33 simulations. The normalized wave number K = 2k x , is displayed as a function of the normalized frequency,W = ω�t , in Fig. 2.
As shown in Fig. 2, the dispersion relation of our LBM method is in good consistency with the exact analytical curve. Additionally, Fig. 2 suggests that the new LBM method has lower dispersion error than the FDTD method for high wave numbers ( K > 2.0) , which is further advantage of the LBM over FDTD method for the specific problem being discussed. On the other hand, the 1D LBM method uses three distribution functions ( f 0 (x, t) , f 1 (x, t), f 2 (x, t) ) to represent two physical fields (H and E), thus requiring ~ 33% more memory compared to the 1D FDTD method. EM wave in Debye medium based on LBM method. In this case, we considered an electrical pulse propagating from air into water. In this case, water is treated as a Debye medium. For water, the properties for the relative permittivity are well described by the Debye model 31,34 with ε s = 81 , ε ∞ = 1.8 , and t 0 = 9.4 × 10 −12 (Fig. 2). It is worth to be noted that, from the Debye model 31,34 (1)  Equation (4) shows that R[ε(ω)] decreases as the frequency of EM wave increases from low frequency to high frequency, whereas Im[ε(ω)] exhibits a bell shape while the frequency changes from low frequency to high frequency. Furthermore, the peak value of Im[ε(ω)] could be calculated as follows: where ω maxε is the frequency that gives the maximum Im[ε(ω)] . For water, as shown in Fig. 3, t 0 = 9.4 × 10 −12 , which gives ω maxε =0.1063 THz. Figure 4 depicts the initial electric field in the air. The computational domain consists of 800 lattices, with the left 400 for free space and the right 400 occupied by water. We used an incident pulse with a maximum amplitude of 1 kV/m and a width of 30 cells. The lattice size was set to 3.75 µm , and the time step was 0.125 ps. These parameters have been used by Luebbers et al. 18 and Chen et al. 12 In Fig. 5, we compared our numerical results with previous studies generated by the well-established (FD) 2 TD and the LBM with pseudo permittivity schemes. It is notable that our results are more consistent with the accurate (FD) 2 TD results than the pseudo permittivity approach.
When dealing with conservative differential equations, it is important to give an estimate of the order of convergence of the model. This can be realized by monitoring the accuracy of the physical variables, specifically, the EM energy density, and the resolution of the lattice grid is increased while keeping the time resolution constant:  www.nature.com/scientificreports/ Figure 6 depicts the EM energy density U as a function of the number of LBM grid cells. We used Richardson's method 30 to evaluate the convergence error of the LBM method, while the exact solution of the EM energy density is given by: with an n + 1 order of error L(δx n+1 ) , where δx = L N , as L is the length of the computational domain, and N is the number of cells. Therefore, the error of the model is at the second-order, and the relative errors can be evaluated as Figure 5 shows that the relative error decreases with an increase in LBM grid cell number as δx 1.92 . The decrease in error, supports that the present scheme has a second-order convergence.

Effects of film thickness and frequency on Gaussian modulated sinusoidal EM pulse. In light
of the validations results presented in sect. "EM wave in Debye medium based on LBM method", the adaptability of the LBM-SEF model is further assessed by examining two intriguing phenomena of the response of dispersive medium to an incoming EM wave. Here we focus on a scenario where an EM pulse travels through the air and encounters a Debye medium. As depicted in Fig. 6, the simulation domain consists of 10,000 lattice cells, each  Here E M is the electric field amplitude of the pulse,ω is 0.2 or 0.8 THz for the low-and high-frequency cases, respectively. The pulse incidents from the air (right), pass through the Debye film (middle), and return back to air, as depicted in Fig. 7.
Depending on the file thickness, the EM wave can be partially or completely aborted by the Debye medium as it passes through it. The simulation results for thin-and thick-films of Debye media are presented in Fig. 8. Figure 8a1,a2 and b1,b2 show that, at time = 0 ps, the EM wave is in the air domain ( Fig. 8a1 and b1). Then, it touches the air-water interface at approximately time = 33 ps ( Fig. 8a2 and b2). At time = 56 ps, Fig. 8a3 shows that the EM wave is partially absorbed by the thin water film, while Fig. 8b3 shows that the EM wave is almost completely absorbed by the thick water film.
The frequency also significantly affects the EM wave propagation. The distribution of the electric fields of the low-and high-frequency Gaussian modulated sinusoidal pulses are depicted in Fig. 9. Figure 9a1 and b1 show that the EM wave at low frequency is in air domain at time = 0 ps, and it then starts touching the air-water interface at time = 33 ps (Fig. 9a2,b2). At time t = 56 ps, the EM wave is split into two constituents: the left in the air is reflected by the Debye medium and travels to the left, while the transmitted wave travels to the right. A comparison between Fig. 9a3 and b3 indicates that the Debye film is more transparent for higher-frequency EM waves. The decrease of the transparency with the EM frequency is relatively low. These results further validate the LBM-SEF method for EM wave propagation in dispersive media.

Conclusions
In this manuscript, we proposed and implemented a novel lattice Boltzmann method with an extended force term to simulate the EM wave propagation in dispersive media. The proposed method is mathematically consistent with Maxwell's equations and the accuracies were demonstrated by comparing simulation results with those of   www.nature.com/scientificreports/ the (FD) 2 TD method and the previous lattice Boltzmann method using pseudo permittivity. Two typical yet intriguing cases were employed to assess the ability and adaptability of the new method for EM waves traveling through Debye films. The results confirmed the accuracy and robustness of the new method for dealing with frequency-dependent reflection and transmission at the air-water interface. Furthermore, the results can aid in constructing two-dimensional and three-dimensional lattice Boltzmann method for the simulation of EM waves in dispersive media.

Materials and methods
The Lattice Boltzmann method (LBM) is widely used to describe the evolution of particle distributions with discrete space and time coordiantes 19 . The density distribution and the particle velocities can be well reproduced by the density and moments of the model 20,21 . Here we briefly describe the technical details of the extended LBM equation with a force term to describe EM wave propagation in dispersive media. We also demonstrated that the new method is consistent with Maxwell equations by using the Chapman-Enskog expansion technique.
Governing equations. We start from the time-dependent Maxwell equations in dielectric media at position x and time t 28 : where H and E denote the intensities of the magnetic and electric fields, respectively, B and D are the magnetic induction and electric displacement,and J is the current density. B and H , D and E are related by B = µH , and D = εE, where µ and ε are the permittivity and permeability of the medium, respectively. Furthermore, µ and ε can also be related to the relative constants µ r ,ε r by µ=µ r µ 0 and ε=ε r ε 0 . As discussed in the study by Chen et al. 12 , the permittivity of a linear, isotropic dispersive medium, ε(ω) , relates D and E in the frequency domain via D(ω) = ε(ω)E(ω) . While in the time domain, we have Here ε ∞ is the permittivity at infinite high frequency, ζ(τ ) is the time domain electric susceptibility. Different from the work by Chen et al. 12 , in this study, we use the following expression: In the previous method by Chen et al. 12 , it is necessary to evaluate the pseudo permittivity for the simulation of the EM wave propagation, which significantly influences the overall accuracy of the simulation via the selection of the time resolution t 12 . However, according to Eq. (12), the evaluation of the pseudo permittivity is no longer required in our method. As will be demonstrated in the following, only a single force term would be sufficient for the simulation.
In a one dimensional dispersive medium, the model can be further simplified for plane EM wave propagation: It is noted that with Eq. (4), the simplified LBM model is also mathematically more strict.
Extended lattice Boltzmann model. In this section, we proposed a dimensionless lattice Boltzmann model with a force term for EM waves in a linear dispersive medium. The proposed lattice Boltzmann model is as follows: where f i (x, t) is the particle distribution function at position x, time t in the i-th direction. e i is the lattice vector in the i-th direction, t is the time step, τ is the relaxation time, and b indexes the lattice vectors. Obviously, b = 2 for one-dimensional media. f eq i (x, t) is the local distribution at equilibrium, and F i (x, t) is a force term. The macroscopic EM characteristics are given by: www.nature.com/scientificreports/ To solve the EM wave equation, we employed the LBM scheme to calculate the bandgap of photonic materials. Furthermore, the f eq i (x, t) is written as: where A i and B i are weighting parameters. According to symmetrical properties, the distribution function at equilibrium, f eq i (x, t) can be written as: Another two conservation constraints were used to specify the distribution weights: Based on Eqs. (17) and (18), we obtain the following equations: Chapman-Enskog expansion. The distribution function f i (x, t) can be expanded up to the second-order according to the Chapman-Enskog expansion: where the expansion parameter θ is usually understood as Knudsen number [40], while it is not necessary to calculate this parameter using the Chapman-Enskog expansion technique.
It shouble be noted that in Eq. (20), f 1 i (x, t) and f 2 i (x, t) can be viewed as formal expansion functions corresponding to distribution functions at different scales. Combining Eqs. (18), (19) and (20), we obtain the following equations: Given that e 0 =0and θ is an infinitesimal, we arrive at Equation (14) can be Taylor expanded with respect to time, The force term in Eq. (14) has the following form: Using Chapman-Enskog expansion, we equate the same orders of ∂ t and ∂ xα , which gives : Assuming the force term takes the form of:  (25) and (26), and substitute into Eq. (14), we have: Grouping terms with same orders of θ and θ 2 : Summing Eq. (28) over i, and using Eqs. (15) and (22), we could get that: and By multipling Eqs. (29) and (30) by e iα and considering (22) we obtain: Comparing the above equation with the original Maxwell Eq. (4), we have the constraints: By considering (θ × Eq. (29) + θ 2 × Eq. (30)), we obtain: Based on Eq. (34), it is evident that the proposed LBM scheme with the force term is capable of recovering the Maxwell equations at the limits of t and θ approach to zero. Furthermore, the force term F i (1) (x, t) can be evaluated by the integration where the force kernel function P(x, t) = t 0 ∂ t E(x, t − τ)ζ(τ)dτ , which will be calculated in sect. "Calculation of the force term of Debye media". is given by: where ε s denotes the static permittivity, ε ∞ denotes the infinite frequency permittivity, t 0 denotes the relaxation time constant, and S(t) denotes the unit step function.
Based on Eqs. (27) and (28), the forcing term is determined by P(x, t) as follows: Furthermore, the term P(x, t) can be calculated based on the following.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.